Packages
if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse") #tidyr, ggplot, dplyr, etc
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr")
if (!require("lme4")) install.packages("lme4")
library("lme4") # mixed-effect models
if (!require("lmerTest")) install.packages("lmerTest")
library("lmerTest") # p-values
if (!require("car")) install.packages("car")
library("car") # test for VIF
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # model export
if (!require("emmeans")) install.packages('emmeans')
library('emmeans')
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # model export
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown formatBackground
This code walks through the statistical models and figures for a study on the rate of change of cutaneous evaporative water loss (CEWL) in response to temperature change in Western Fence Lizards (Sceloporus occidentalis) (published in the Journal __ in 2023). Data was collected and analyzed by Calvin Davis and Savannah Weaver, under the advising of Dr. Emily Taylor at California Polytechnic State University.
Load Data
Get the RDS files created in the data wrangling Rmd:
temp_time_series <- read_rds("Data/temp_time_series.RDS")
CEWL_time_series <- read_rds("Data/CEWL_time_series.RDS")Stats & Models
Temp Correlation
Check the correlation between dorsal and cloacal temps:
temp_corr <- lm(data = temp_time_series,
calibrated_dorsal_temp ~ calibrated_cloacal_temp)
summary(temp_corr)##
## Call:
## lm(formula = calibrated_dorsal_temp ~ calibrated_cloacal_temp,
## data = temp_time_series)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0389 -0.7784 0.0075 0.5911 5.3351
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.966310 0.049820 139.8 <2e-16 ***
## calibrated_cloacal_temp 0.747869 0.001899 393.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.388 on 25249 degrees of freedom
## Multiple R-squared: 0.86, Adjusted R-squared: 0.8599
## F-statistic: 1.55e+05 on 1 and 25249 DF, p-value: < 2.2e-16#write.csv(broom.mixed::tidy(temp_corr),
# "./Results_Statistics/body_temp_correlation.csv")They are actually pretty different, so run parallel models of CEWL, one with dorsal temp, and one with cloacal temp.
CEWL (raw) ~ temp LMM
First, test linear effects of temperature:
CEWL_LMM_02a <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
calibrated_cloacal_temp *
treatment +
(1|date/lizard_ID))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluessummary(CEWL_LMM_02a)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 87927.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1047 -0.5895 -0.0371 0.4647 8.7021
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 10.066 3.173
## date (Intercept) 22.540 4.748
## Residual 1.882 1.372
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 51.20985 2.56934 19.93
## calibrated_cloacal_temp -1.31478 0.03945 -33.33
## treatmentCooling -51.86569 1.71910 -30.17
## treatmentHeating -39.66638 1.78910 -22.17
## calibrated_cloacal_temp:treatmentCooling 1.62627 0.03957 41.09
## calibrated_cloacal_temp:treatmentHeating 1.54114 0.03957 38.95
##
## Correlation of Fixed Effects:
## (Intr) clbr__ trtmnC trtmnH cl__:C
## clbrtd_clc_ -0.421
## tretmntClng -0.472 0.628
## tretmntHtng -0.459 0.605 0.669
## clbrtd_c_:C 0.419 -0.997 -0.629 -0.603
## clbrtd_c_:H 0.419 -0.997 -0.626 -0.607 0.994
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesvif(CEWL_LMM_02a)## GVIF Df GVIF^(1/(2*Df))
## calibrated_cloacal_temp 309.047830 1 17.579756
## treatment 1.853502 2 1.166805
## calibrated_cloacal_temp:treatment 311.665946 2 4.201674CEWL_LMM_02b <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
calibrated_dorsal_temp *
treatment +
(1|date/lizard_ID))
summary(CEWL_LMM_02b)## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 88772.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9809 -0.5887 -0.0570 0.4961 8.9013
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 10.316 3.212
## date (Intercept) 21.015 4.584
## Residual 1.946 1.395
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 51.27388 2.54285 20.16
## calibrated_dorsal_temp -1.30362 0.04142 -31.48
## treatmentCooling -53.18221 1.77431 -29.97
## treatmentHeating -40.88244 1.84485 -22.16
## calibrated_dorsal_temp:treatmentCooling 1.64631 0.04157 39.60
## calibrated_dorsal_temp:treatmentHeating 1.57709 0.04164 37.88
##
## Correlation of Fixed Effects:
## (Intr) clbr__ trtmnC trtmnH cl__:C
## clbrtd_drs_ -0.451
## tretmntClng -0.499 0.645
## tretmntHtng -0.486 0.623 0.680
## clbrtd_d_:C 0.449 -0.996 -0.647 -0.620
## clbrtd_d_:H 0.449 -0.995 -0.642 -0.625 0.991vif(CEWL_LMM_02b)## GVIF Df GVIF^(1/(2*Df))
## calibrated_dorsal_temp 231.169627 1 15.204263
## treatment 1.944891 2 1.180929
## calibrated_dorsal_temp:treatment 234.207111 2 3.912011These are REALLY good models, but based on the figures below, the relationship isn’t actually linear… so test the inclusion of polynomial factors:
CEWL_LMM_03a <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
poly(calibrated_cloacal_temp, 2)*treatment +
(1|date/lizard_ID))
summary(CEWL_LMM_03a)## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 81611.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6892 -0.6401 -0.0881 0.5378 7.8576
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 11.773 3.431
## date (Intercept) 23.794 4.878
## Residual 1.468 1.212
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 0.6557 2.4503
## poly(calibrated_cloacal_temp, 2)1 1422.2807 62.9113
## poly(calibrated_cloacal_temp, 2)2 -3419.1801 82.4963
## treatmentCooling 6.6956 1.5010
## treatmentHeating 16.7664 1.5897
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling -1162.0238 62.9471
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling 3565.6587 82.5285
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating -1223.6577 62.9510
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating 3338.2055 82.5256
## t value
## (Intercept) 0.268
## poly(calibrated_cloacal_temp, 2)1 22.608
## poly(calibrated_cloacal_temp, 2)2 -41.446
## treatmentCooling 4.461
## treatmentHeating 10.547
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling -18.460
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling 43.205
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating -19.438
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating 40.451
##
## Correlation of Fixed Effects:
## (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.159
## ply(c__,2)2 0.164 -0.914
## tretmntClng -0.336 0.260 -0.268
## tretmntHtng -0.324 0.244 -0.251 0.505
## pl(__,2)1:C 0.159 -0.999 0.914 -0.259 -0.244
## pl(__,2)2:C -0.164 0.914 -1.000 0.268 0.251 -0.913
## pl(__,2)1:H 0.159 -0.999 0.914 -0.259 -0.244 0.999 -0.913
## pl(__,2)2:H -0.164 0.914 -1.000 0.268 0.251 -0.914 0.999
## p(__,2)1:H
## ply(c__,2)1
## ply(c__,2)2
## tretmntClng
## tretmntHtng
## pl(__,2)1:C
## pl(__,2)2:C
## pl(__,2)1:H
## pl(__,2)2:H -0.914vif(CEWL_LMM_03a)## GVIF Df GVIF^(1/(2*Df))
## poly(calibrated_cloacal_temp, 2) 9.312576e+05 2 31.064721
## treatment 1.100911e+00 2 1.024326
## poly(calibrated_cloacal_temp, 2):treatment 9.312225e+05 4 5.573547anova(CEWL_LMM_03a, CEWL_LMM_02a)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_02a: CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_02a 9 87929 88002 -43956 87911
## CEWL_LMM_03a 12 81676 81774 -40826 81652 6259.3 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1CEWL_LMM_03b <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
poly(calibrated_dorsal_temp, 2)*treatment +
(1|date/lizard_ID))
summary(CEWL_LMM_03b)## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 83028.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7628 -0.6345 -0.0833 0.5160 8.3769
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 9.996 3.162
## date (Intercept) 21.953 4.685
## Residual 1.554 1.246
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.248 2.336 1.390
## poly(calibrated_dorsal_temp, 2)1 1534.609 67.430 22.759
## poly(calibrated_dorsal_temp, 2)2 -2711.185 75.106 -36.098
## treatmentCooling 3.529 1.386 2.546
## treatmentHeating 14.341 1.468 9.766
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling -1350.765 67.457 -20.024
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling 2852.577 75.134 37.967
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating -1365.961 67.469 -20.246
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating 2666.790 75.141 35.490
##
## Correlation of Fixed Effects:
## (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.162
## ply(c__,2)2 0.163 -0.946
## tretmntClng -0.327 0.273 -0.276
## tretmntHtng -0.315 0.257 -0.260 0.507
## pl(__,2)1:C 0.162 -1.000 0.946 -0.273 -0.257
## pl(__,2)2:C -0.163 0.946 -1.000 0.276 0.259 -0.946
## pl(__,2)1:H 0.162 -0.999 0.946 -0.273 -0.257 0.999 -0.945
## pl(__,2)2:H -0.163 0.946 -1.000 0.276 0.259 -0.945 0.999
## p(__,2)1:H
## ply(c__,2)1
## ply(c__,2)2
## tretmntClng
## tretmntHtng
## pl(__,2)1:C
## pl(__,2)2:C
## pl(__,2)1:H
## pl(__,2)2:H -0.945vif(CEWL_LMM_03b)## GVIF Df GVIF^(1/(2*Df))
## poly(calibrated_dorsal_temp, 2) 5.528187e+05 2 27.267523
## treatment 1.107683e+00 2 1.025897
## poly(calibrated_dorsal_temp, 2):treatment 5.527949e+05 4 5.221803anova(CEWL_LMM_03b, CEWL_LMM_02b)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_02b: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_02b 9 88775 88848 -44379 88757
## CEWL_LMM_03b 12 83092 83190 -41534 83068 5689.1 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The poly models have significantly better explanatory power than the linear models.
Also try log transformation to make sure a polynomial effect is best:
CEWL_LMM_04a <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
(log(calibrated_cloacal_temp))*treatment +
(1|date/lizard_ID))
summary(CEWL_LMM_04a)## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 88908.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1896 -0.5896 -0.0624 0.4489 8.9787
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 9.877 3.143
## date (Intercept) 22.480 4.741
## Residual 1.958 1.399
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 129.257 4.312 29.97
## log(calibrated_cloacal_temp) -34.468 1.097 -31.41
## treatmentCooling -143.864 3.871 -37.16
## treatmentHeating -131.365 3.906 -33.63
## log(calibrated_cloacal_temp):treatmentCooling 41.258 1.100 37.50
## log(calibrated_cloacal_temp):treatmentHeating 40.503 1.101 36.80
##
## Correlation of Fixed Effects:
## (Intr) lg(__) trtmnC trtmnH l(__):C
## lg(clbrt__) -0.842
## tretmntClng -0.844 0.938
## tretmntHtng -0.838 0.930 0.930
## lg(clb__):C 0.840 -0.998 -0.940 -0.928
## lg(clb__):H 0.840 -0.997 -0.935 -0.933 0.995vif(CEWL_LMM_04a)## GVIF Df GVIF^(1/(2*Df))
## log(calibrated_cloacal_temp) 382.13042 1 19.548156
## treatment 11.39233 2 1.837186
## log(calibrated_cloacal_temp):treatment 416.69458 2 4.518086anova(CEWL_LMM_04a, CEWL_LMM_02a)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_04a: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_02a: CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04a 9 88930 89003 -44456 88912
## CEWL_LMM_02a 9 87929 88002 -43956 87911 1000.4 0anova(CEWL_LMM_04a, CEWL_LMM_03a)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_04a: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04a 9 88930 89003 -44456 88912
## CEWL_LMM_03a 12 81676 81774 -40826 81652 7259.7 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1CEWL_LMM_04b <- lme4::lmer(data = temp_time_series,
CEWL_g_m2h ~
log(calibrated_dorsal_temp)*treatment +
treatment +
(1|date/lizard_ID))
summary(CEWL_LMM_04b)## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment +
## (1 | date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 89739.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0661 -0.5838 -0.0705 0.4871 9.0055
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 10.302 3.210
## date (Intercept) 21.092 4.593
## Residual 2.024 1.423
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 131.771 4.497 29.30
## log(calibrated_dorsal_temp) -35.114 1.169 -30.05
## treatmentCooling -151.142 4.117 -36.71
## treatmentHeating -137.198 4.155 -33.02
## log(calibrated_dorsal_temp):treatmentCooling 43.244 1.172 36.89
## log(calibrated_dorsal_temp):treatmentHeating 42.178 1.174 35.94
##
## Correlation of Fixed Effects:
## (Intr) lg(__) trtmnC trtmnH l(__):C
## lg(clbrt__) -0.863
## tretmntClng -0.863 0.942
## tretmntHtng -0.857 0.934 0.933
## lg(clb__):C 0.860 -0.997 -0.945 -0.931
## lg(clb__):H 0.859 -0.996 -0.938 -0.938 0.993vif(CEWL_LMM_04b)## GVIF Df GVIF^(1/(2*Df))
## log(calibrated_dorsal_temp) 275.88918 1 16.609912
## treatment 12.81722 2 1.892119
## log(calibrated_dorsal_temp):treatment 315.47270 2 4.214446anova(CEWL_LMM_04b, CEWL_LMM_02b)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_04b: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_02b: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04b 9 89762 89835 -44872 89744
## CEWL_LMM_02b 9 88775 88848 -44379 88757 987.18 0anova(CEWL_LMM_04b, CEWL_LMM_03b)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## CEWL_LMM_04b: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04b 9 89762 89835 -44872 89744
## CEWL_LMM_03b 12 83092 83190 -41534 83068 6676.3 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The logarithmic models are NOT better than the linear or polynomial ones.
VIFs are still close to 1 with the polynomial included, and all variables seem significant, so check assumptions:
plot(CEWL_LMM_03a)
plot(CEWL_LMM_03b)
They both have a really weird spike, but otherwise there’s no fanning or different pattern.
Re-run with p-values:
CEWL_LMM_besta <- lmerTest::lmer(data = temp_time_series,
CEWL_g_m2h ~
poly(calibrated_cloacal_temp, 2) * treatment +
(1|date/lizard_ID))
summary(CEWL_LMM_besta)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 81611.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6892 -0.6401 -0.0881 0.5378 7.8576
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 11.773 3.431
## date (Intercept) 23.794 4.878
## Residual 1.468 1.212
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 0.6557 2.4503
## poly(calibrated_cloacal_temp, 2)1 1422.2807 62.9113
## poly(calibrated_cloacal_temp, 2)2 -3419.1801 82.4963
## treatmentCooling 6.6956 1.5010
## treatmentHeating 16.7664 1.5897
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling -1162.0238 62.9471
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling 3565.6587 82.5285
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating -1223.6577 62.9510
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating 3338.2055 82.5256
## df t value Pr(>|t|)
## (Intercept) 5.4550 0.268 0.798840
## poly(calibrated_cloacal_temp, 2)1 25235.9702 22.608 < 2e-16
## poly(calibrated_cloacal_temp, 2)2 25238.1540 -41.446 < 2e-16
## treatmentCooling 30.3218 4.461 0.000104
## treatmentHeating 29.8965 10.547 1.36e-11
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling 25235.9622 -18.460 < 2e-16
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling 25238.1521 43.205 < 2e-16
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating 25235.9530 -19.438 < 2e-16
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating 25238.1526 40.451 < 2e-16
##
## (Intercept)
## poly(calibrated_cloacal_temp, 2)1 ***
## poly(calibrated_cloacal_temp, 2)2 ***
## treatmentCooling ***
## treatmentHeating ***
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling ***
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling ***
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating ***
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.159
## ply(c__,2)2 0.164 -0.914
## tretmntClng -0.336 0.260 -0.268
## tretmntHtng -0.324 0.244 -0.251 0.505
## pl(__,2)1:C 0.159 -0.999 0.914 -0.259 -0.244
## pl(__,2)2:C -0.164 0.914 -1.000 0.268 0.251 -0.913
## pl(__,2)1:H 0.159 -0.999 0.914 -0.259 -0.244 0.999 -0.913
## pl(__,2)2:H -0.164 0.914 -1.000 0.268 0.251 -0.914 0.999
## p(__,2)1:H
## ply(c__,2)1
## ply(c__,2)2
## tretmntClng
## tretmntHtng
## pl(__,2)1:C
## pl(__,2)2:C
## pl(__,2)1:H
## pl(__,2)2:H -0.914#anova(CEWL_LMM_besta, type = "3", ddf = "Kenward-Roger")
CEWL_LMM_bestb <- lmerTest::lmer(data = temp_time_series,
CEWL_g_m2h ~
poly(calibrated_dorsal_temp, 2) * treatment +
(1|date/lizard_ID))
summary(CEWL_LMM_bestb)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 83028.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7628 -0.6345 -0.0833 0.5160 8.3769
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 9.996 3.162
## date (Intercept) 21.953 4.685
## Residual 1.554 1.246
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 3.248 2.336
## poly(calibrated_dorsal_temp, 2)1 1534.609 67.430
## poly(calibrated_dorsal_temp, 2)2 -2711.185 75.106
## treatmentCooling 3.529 1.386
## treatmentHeating 14.341 1.468
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling -1350.765 67.457
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling 2852.577 75.134
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating -1365.961 67.469
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating 2666.790 75.141
## df t value Pr(>|t|)
## (Intercept) 5.353 1.390 0.2194
## poly(calibrated_dorsal_temp, 2)1 25237.527 22.759 < 2e-16
## poly(calibrated_dorsal_temp, 2)2 25236.097 -36.098 < 2e-16
## treatmentCooling 30.609 2.546 0.0162
## treatmentHeating 30.174 9.766 7.46e-11
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling 25237.523 -20.024 < 2e-16
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling 25236.096 37.967 < 2e-16
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating 25237.517 -20.246 < 2e-16
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating 25236.103 35.490 < 2e-16
##
## (Intercept)
## poly(calibrated_dorsal_temp, 2)1 ***
## poly(calibrated_dorsal_temp, 2)2 ***
## treatmentCooling *
## treatmentHeating ***
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling ***
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling ***
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating ***
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.162
## ply(c__,2)2 0.163 -0.946
## tretmntClng -0.327 0.273 -0.276
## tretmntHtng -0.315 0.257 -0.260 0.507
## pl(__,2)1:C 0.162 -1.000 0.946 -0.273 -0.257
## pl(__,2)2:C -0.163 0.946 -1.000 0.276 0.259 -0.946
## pl(__,2)1:H 0.162 -0.999 0.946 -0.273 -0.257 0.999 -0.945
## pl(__,2)2:H -0.163 0.946 -1.000 0.276 0.259 -0.945 0.999
## p(__,2)1:H
## ply(c__,2)1
## ply(c__,2)2
## tretmntClng
## tretmntHtng
## pl(__,2)1:C
## pl(__,2)2:C
## pl(__,2)1:H
## pl(__,2)2:H -0.945#anova(CEWL_LMM_bestb, type = "3", ddf = "Kenward-Roger")
# double check sample sizes
temp_time_series %>%
group_by(lizard_ID, treatment) %>%
summarise(n()) %>%
group_by(treatment) %>%
summarise(n())## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.## # A tibble: 3 × 2
## treatment `n()`
## <fct> <int>
## 1 Control 11
## 2 Cooling 12
## 3 Heating 10These are the best models for CEWL ~ temperature.
broom.mixed::tidy(CEWL_LMM_besta) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/CEWL_best_model_clotemp.csv")
broom.mixed::tidy(CEWL_LMM_bestb) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/CEWL_best_model_dorstemp.csv")CEWL (raw) ~ time LMM
NOTE: Start times vary from 97-170 seconds after starting time series.
Make a polynomial model first:
rates_LMM_01 <- lme4::lmer(data = CEWL_time_series,
CEWL_g_m2h ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
summary(rates_LMM_01)## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## Data: CEWL_time_series
##
## REML criterion at convergence: 96853.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3325 -0.5700 -0.0438 0.5178 8.0025
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 9.976 3.158
## date (Intercept) 17.817 4.221
## Residual 1.754 1.324
## Number of obs: 28405, groups: lizard_ID:date, 37; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.957 2.099 7.126
## poly(time_min, 2)1 -137.974 2.309 -59.765
## poly(time_min, 2)2 51.100 2.322 22.008
## treatmentCooling -8.473 1.297 -6.533
## treatmentHeating 3.416 1.280 2.669
## poly(time_min, 2)1:treatmentCooling -88.851 3.376 -26.319
## poly(time_min, 2)2:treatmentCooling 83.367 3.337 24.980
## poly(time_min, 2)1:treatmentHeating 301.172 3.201 94.092
## poly(time_min, 2)2:treatmentHeating -108.483 3.201 -33.893
##
## Correlation of Fixed Effects:
## (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 0.000
## ply(tm_,2)2 0.000 -0.067
## tretmntClng -0.305 0.001 0.000
## tretmntHtng -0.317 0.001 0.000 0.499
## ply(_,2)1:C 0.000 -0.684 0.046 0.000 0.000
## ply(_,2)2:C 0.000 0.047 -0.696 0.000 0.000 0.008
## ply(_,2)1:H 0.000 -0.721 0.049 0.000 0.000 0.493 -0.034
## ply(_,2)2:H 0.000 0.049 -0.725 0.000 0.000 -0.033 0.505
## p(_,2)1:H
## ply(tm_,2)1
## ply(tm_,2)2
## tretmntClng
## tretmntHtng
## ply(_,2)1:C
## ply(_,2)2:C
## ply(_,2)1:H
## ply(_,2)2:H -0.029vif(rates_LMM_01)## GVIF Df GVIF^(1/(2*Df))
## poly(time_min, 2) 9.050897 2 1.734494
## treatment 1.000002 2 1.000000
## poly(time_min, 2):treatment 9.050904 4 1.317002drop1(rates_LMM_01)## Single term deletions
##
## Model:
## CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC
## <none> 96906
## poly(time_min, 2):treatment 4 111987Also make a linear model version:
rates_LMM_02 <- lme4::lmer(data = CEWL_time_series,
CEWL_g_m2h ~
time_min * treatment +
(1|date/lizard_ID))
summary(rates_LMM_02)## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
## Data: CEWL_time_series
##
## REML criterion at convergence: 100915.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6381 -0.5611 -0.0446 0.5129 8.5948
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 10.08 3.175
## date (Intercept) 17.74 4.212
## Residual 2.02 1.421
## Number of obs: 28405, groups: lizard_ID:date, 37; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 16.620731 2.097842 7.923
## time_min -0.205350 0.003773 -54.431
## treatmentCooling -7.212880 1.304581 -5.529
## treatmentHeating -0.277100 1.287525 -0.215
## time_min:treatmentCooling -0.156727 0.005515 -28.418
## time_min:treatmentHeating 0.455474 0.005236 86.986
##
## Correlation of Fixed Effects:
## (Intr) tim_mn trtmnC trtmnH tm_m:C
## time_min -0.015
## tretmntClng -0.307 0.024
## tretmntHtng -0.319 0.024 0.500
## tm_mn:trtmC 0.010 -0.684 -0.034 -0.017
## tm_mn:trtmH 0.011 -0.720 -0.017 -0.033 0.493anova(rates_LMM_02, rates_LMM_01)## refitting model(s) with ML (instead of REML)## Data: CEWL_time_series
## Models:
## rates_LMM_02: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rates_LMM_02 9 100913 100987 -50447 100895
## rates_LMM_01 12 96906 97005 -48441 96882 4012.3 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Polynomial is much better than linear.
Also compare to log version:
rates_LMM_03 <- lme4::lmer(data = CEWL_time_series,
CEWL_g_m2h ~
log(time_min) * treatment +
(1|date/lizard_ID))
summary(rates_LMM_03)## Linear mixed model fit by REML ['lmerMod']
## Formula: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
## Data: CEWL_time_series
##
## REML criterion at convergence: 97113.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1092 -0.5865 -0.0434 0.5019 7.9036
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 10.057 3.171
## date (Intercept) 17.764 4.215
## Residual 1.767 1.329
## Number of obs: 28405, groups: lizard_ID:date, 37; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 17.83671 2.09858 8.499
## log(time_min) -1.48000 0.02353 -62.902
## treatmentCooling -6.24930 1.30381 -4.793
## treatmentHeating -2.77380 1.28675 -2.156
## log(time_min):treatmentCooling -1.15769 0.03356 -34.500
## log(time_min):treatmentHeating 3.18416 0.03220 98.902
##
## Correlation of Fixed Effects:
## (Intr) lg(t_) trtmnC trtmnH l(_):C
## log(tim_mn) -0.022
## tretmntClng -0.307 0.036
## tretmntHtng -0.319 0.036 0.500
## lg(tm_mn):C 0.016 -0.701 -0.050 -0.025
## lg(tm_mn):H 0.016 -0.731 -0.026 -0.049 0.512anova(rates_LMM_02, rates_LMM_03)## refitting model(s) with ML (instead of REML)## Data: CEWL_time_series
## Models:
## rates_LMM_02: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
## rates_LMM_03: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rates_LMM_02 9 100913 100987 -50447 100895
## rates_LMM_03 9 97121 97196 -48552 97103 3791.4 0anova(rates_LMM_01, rates_LMM_03)## refitting model(s) with ML (instead of REML)## Data: CEWL_time_series
## Models:
## rates_LMM_03: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rates_LMM_03 9 97121 97196 -48552 97103
## rates_LMM_01 12 96906 97005 -48441 96882 220.96 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The polynomial model is best (better than both log and lm options).
Check linear regression assumptions:
plot(rates_LMM_01)
Again, there is a really weird spike, but otherwise there’s no fanning or pattern.
re-run with p-values:
rates_LMM_best <- lmerTest::lmer(data = CEWL_time_series,
CEWL_g_m2h ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
summary(rates_LMM_best)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## Data: CEWL_time_series
##
## REML criterion at convergence: 96853.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3325 -0.5700 -0.0438 0.5178 8.0025
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 9.976 3.158
## date (Intercept) 17.817 4.221
## Residual 1.754 1.324
## Number of obs: 28405, groups: lizard_ID:date, 37; date, 5
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 14.957 2.099 5.189 7.126
## poly(time_min, 2)1 -137.974 2.309 28362.024 -59.765
## poly(time_min, 2)2 51.100 2.322 28362.016 22.008
## treatmentCooling -8.473 1.297 30.005 -6.533
## treatmentHeating 3.416 1.280 30.061 2.669
## poly(time_min, 2)1:treatmentCooling -88.851 3.376 28362.571 -26.319
## poly(time_min, 2)2:treatmentCooling 83.367 3.337 28362.083 24.980
## poly(time_min, 2)1:treatmentHeating 301.172 3.201 28362.093 94.092
## poly(time_min, 2)2:treatmentHeating -108.483 3.201 28362.048 -33.893
## Pr(>|t|)
## (Intercept) 0.000723 ***
## poly(time_min, 2)1 < 2e-16 ***
## poly(time_min, 2)2 < 2e-16 ***
## treatmentCooling 3.17e-07 ***
## treatmentHeating 0.012151 *
## poly(time_min, 2)1:treatmentCooling < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 0.000
## ply(tm_,2)2 0.000 -0.067
## tretmntClng -0.305 0.001 0.000
## tretmntHtng -0.317 0.001 0.000 0.499
## ply(_,2)1:C 0.000 -0.684 0.046 0.000 0.000
## ply(_,2)2:C 0.000 0.047 -0.696 0.000 0.000 0.008
## ply(_,2)1:H 0.000 -0.721 0.049 0.000 0.000 0.493 -0.034
## ply(_,2)2:H 0.000 0.049 -0.725 0.000 0.000 -0.033 0.505
## p(_,2)1:H
## ply(tm_,2)1
## ply(tm_,2)2
## tretmntClng
## tretmntHtng
## ply(_,2)1:C
## ply(_,2)2:C
## ply(_,2)1:H
## ply(_,2)2:H -0.029#anova(rates_LMM_best, type = "3", ddf = "Kenward-Roger")
# double check sample sizes
CEWL_time_series %>%
group_by(lizard_ID, treatment) %>%
summarise(n()) %>%
group_by(treatment) %>%
summarise(n())## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.## # A tibble: 3 × 2
## treatment `n()`
## <fct> <int>
## 1 Control 12
## 2 Cooling 12
## 3 Heating 13Treatment explained some variation, and the effect of time explained most of the variation (F(,) = , p <)
broom.mixed::tidy(rates_LMM_best) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/CEWL_best_model_time.csv")Body Temperature ~ Time
linear models:
body_temp_LMM_01a <- lme4::lmer(data = temp_time_series,
calibrated_cloacal_temp ~
time_min * treatment +
(1|date/lizard_ID))
summary(body_temp_LMM_01a)## Linear mixed model fit by REML ['lmerMod']
## Formula: calibrated_cloacal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 63649.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9785 -0.4675 -0.0268 0.3484 6.5452
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 3.0775 1.7543
## date (Intercept) 0.1705 0.4129
## Residual 0.7198 0.8484
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 26.728519 0.561393 47.611
## time_min 0.084958 0.002356 36.058
## treatmentCooling 6.075852 0.734617 8.271
## treatmentHeating -10.116980 0.772484 -13.097
## time_min:treatmentCooling -1.255658 0.003364 -373.236
## time_min:treatmentHeating 1.161071 0.003411 340.363
##
## Correlation of Fixed Effects:
## (Intr) tim_mn trtmnC trtmnH tm_m:C
## time_min -0.035
## tretmntClng -0.680 0.027
## tretmntHtng -0.651 0.025 0.491
## tm_mn:trtmC 0.024 -0.700 -0.037 -0.018
## tm_mn:trtmH 0.024 -0.691 -0.018 -0.036 0.484vif(body_temp_LMM_01a)## GVIF Df GVIF^(1/(2*Df))
## time_min 2.875009 1 1.695585
## treatment 1.002620 2 1.000654
## time_min:treatment 2.879911 2 1.302701body_temp_LMM_01b <- lme4::lmer(data = temp_time_series,
calibrated_dorsal_temp ~
time_min * treatment +
(1|date/lizard_ID))
summary(body_temp_LMM_01b)## Linear mixed model fit by REML ['lmerMod']
## Formula: calibrated_dorsal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 77541.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8874 -0.3920 -0.0119 0.2875 4.8825
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 1.5354 1.2391
## date (Intercept) 0.3511 0.5925
## Residual 1.2497 1.1179
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 27.094037 0.459926 58.91
## time_min 0.077778 0.003105 25.05
## treatmentCooling 6.418579 0.521574 12.31
## treatmentHeating -7.961526 0.551634 -14.43
## time_min:treatmentCooling -1.142811 0.004433 -257.81
## time_min:treatmentHeating 0.786695 0.004495 175.02
##
## Correlation of Fixed Effects:
## (Intr) tim_mn trtmnC trtmnH tm_m:C
## time_min -0.056
## tretmntClng -0.586 0.050
## tretmntHtng -0.562 0.047 0.480
## tm_mn:trtmC 0.039 -0.700 -0.068 -0.033
## tm_mn:trtmH 0.039 -0.691 -0.034 -0.066 0.484vif(body_temp_LMM_01b)## GVIF Df GVIF^(1/(2*Df))
## time_min 2.875083 1 1.695607
## treatment 1.008948 2 1.002230
## time_min:treatment 2.891679 2 1.304030polynomial models:
body_temp_LMM_02a <- lme4::lmer(data = temp_time_series,
calibrated_cloacal_temp ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
summary(body_temp_LMM_02a)## Linear mixed model fit by REML ['lmerMod']
## Formula: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 41254.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8989 -0.4649 -0.0011 0.3942 6.3225
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 3.0570 1.7484
## date (Intercept) 0.1856 0.4308
## Residual 0.2967 0.5447
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 27.4181 0.5621 48.780
## poly(time_min, 2)1 52.0361 0.9376 55.500
## poly(time_min, 2)2 6.0643 0.9435 6.428
## treatmentCooling -4.0686 0.7317 -5.560
## treatmentHeating -0.7262 0.7697 -0.944
## poly(time_min, 2)1:treatmentCooling -764.1781 1.3382 -571.054
## poly(time_min, 2)2:treatmentCooling 161.4913 1.3255 121.838
## poly(time_min, 2)1:treatmentHeating 716.1235 1.3557 528.233
## poly(time_min, 2)2:treatmentHeating -63.6826 1.3531 -47.063
##
## Correlation of Fixed Effects:
## (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001
## ply(tm_,2)2 0.000 -0.077
## tretmntClng -0.676 0.000 0.000
## tretmntHtng -0.647 0.000 0.000 0.490
## ply(_,2)1:C 0.000 -0.701 0.054 0.000 0.000
## ply(_,2)2:C 0.000 0.055 -0.712 0.000 0.000 -0.005
## ply(_,2)1:H 0.000 -0.692 0.053 0.000 0.000 0.485 -0.038
## ply(_,2)2:H 0.000 0.054 -0.697 0.000 0.000 -0.038 0.496
## p(_,2)1:H
## ply(tm_,2)1
## ply(tm_,2)2
## tretmntClng
## tretmntHtng
## ply(_,2)1:C
## ply(_,2)2:C
## ply(_,2)1:H
## ply(_,2)2:H -0.021car::vif(body_temp_LMM_02a)## GVIF Df GVIF^(1/(2*Df))
## poly(time_min, 2) 8.580314 2 1.711496
## treatment 1.000001 2 1.000000
## poly(time_min, 2):treatment 8.580316 4 1.308241anova(body_temp_LMM_02a, body_temp_LMM_01a)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## body_temp_LMM_01a: calibrated_cloacal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## body_temp_LMM_02a: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## body_temp_LMM_01a 9 63638 63712 -31810 63620
## body_temp_LMM_02a 12 41291 41388 -20633 41267 22354 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1body_temp_LMM_02b <- lme4::lmer(data = temp_time_series,
calibrated_dorsal_temp ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
summary(body_temp_LMM_02b)## Linear mixed model fit by REML ['lmerMod']
## Formula: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 71392.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8879 -0.3478 -0.0225 0.3300 6.0970
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 1.5384 1.2403
## date (Intercept) 0.3563 0.5969
## Residual 0.9812 0.9906
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 27.7236 0.4606 60.191
## poly(time_min, 2)1 48.0115 1.7051 28.158
## poly(time_min, 2)2 0.6860 1.7158 0.400
## treatmentCooling -2.8157 0.5208 -5.406
## treatmentHeating -1.5969 0.5509 -2.899
## poly(time_min, 2)1:treatmentCooling -697.0769 2.4336 -286.444
## poly(time_min, 2)2:treatmentCooling 133.9832 2.4104 55.585
## poly(time_min, 2)1:treatmentHeating 484.8482 2.4654 196.660
## poly(time_min, 2)2:treatmentHeating -42.9964 2.4608 -17.473
##
## Correlation of Fixed Effects:
## (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001
## ply(tm_,2)2 0.000 -0.077
## tretmntClng -0.584 0.001 0.000
## tretmntHtng -0.560 0.001 0.000 0.480
## ply(_,2)1:C 0.001 -0.701 0.054 0.000 -0.001
## ply(_,2)2:C 0.000 0.055 -0.712 0.000 0.000 -0.005
## ply(_,2)1:H 0.001 -0.692 0.053 -0.001 0.000 0.485 -0.038
## ply(_,2)2:H 0.000 0.054 -0.697 0.000 0.000 -0.038 0.496
## p(_,2)1:H
## ply(tm_,2)1
## ply(tm_,2)2
## tretmntClng
## tretmntHtng
## ply(_,2)1:C
## ply(_,2)2:C
## ply(_,2)1:H
## ply(_,2)2:H -0.021car::vif(body_temp_LMM_02b)## GVIF Df GVIF^(1/(2*Df))
## poly(time_min, 2) 8.580549 2 1.711507
## treatment 1.000007 2 1.000002
## poly(time_min, 2):treatment 8.580567 4 1.308246anova(body_temp_LMM_02b, body_temp_LMM_01b)## refitting model(s) with ML (instead of REML)## Data: temp_time_series
## Models:
## body_temp_LMM_01b: calibrated_dorsal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## body_temp_LMM_02b: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## body_temp_LMM_01b 9 77531 77604 -38756 77513
## body_temp_LMM_02b 12 71434 71532 -35705 71410 6102.6 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Again, the polynomial models are MUCH better.
body_temp_LMM_besta <- lmerTest::lmer(data = temp_time_series,
calibrated_cloacal_temp ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
summary(body_temp_LMM_besta)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 41254.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8989 -0.4649 -0.0011 0.3942 6.3225
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 3.0570 1.7484
## date (Intercept) 0.1856 0.4308
## Residual 0.2967 0.5447
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 27.4181 0.5621 16.5563 48.780
## poly(time_min, 2)1 52.0361 0.9376 25212.0138 55.500
## poly(time_min, 2)2 6.0643 0.9435 25212.0089 6.428
## treatmentCooling -4.0686 0.7317 26.4065 -5.560
## treatmentHeating -0.7262 0.7697 27.6621 -0.944
## poly(time_min, 2)1:treatmentCooling -764.1781 1.3382 25212.2800 -571.054
## poly(time_min, 2)2:treatmentCooling 161.4913 1.3255 25212.0407 121.838
## poly(time_min, 2)1:treatmentHeating 716.1235 1.3557 25212.0589 528.233
## poly(time_min, 2)2:treatmentHeating -63.6826 1.3531 25212.0326 -47.063
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## poly(time_min, 2)1 < 2e-16 ***
## poly(time_min, 2)2 1.32e-10 ***
## treatmentCooling 7.32e-06 ***
## treatmentHeating 0.354
## poly(time_min, 2)1:treatmentCooling < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001
## ply(tm_,2)2 0.000 -0.077
## tretmntClng -0.676 0.000 0.000
## tretmntHtng -0.647 0.000 0.000 0.490
## ply(_,2)1:C 0.000 -0.701 0.054 0.000 0.000
## ply(_,2)2:C 0.000 0.055 -0.712 0.000 0.000 -0.005
## ply(_,2)1:H 0.000 -0.692 0.053 0.000 0.000 0.485 -0.038
## ply(_,2)2:H 0.000 0.054 -0.697 0.000 0.000 -0.038 0.496
## p(_,2)1:H
## ply(tm_,2)1
## ply(tm_,2)2
## tretmntClng
## tretmntHtng
## ply(_,2)1:C
## ply(_,2)2:C
## ply(_,2)1:H
## ply(_,2)2:H -0.021write.csv(broom.mixed::tidy(body_temp_LMM_besta),
"./Results_Statistics/temp_time_best_model_clotemp.csv")
body_temp_LMM_bestb <- lmerTest::lmer(data = temp_time_series,
calibrated_dorsal_temp ~
poly(time_min, 2) * treatment +
(1|date/lizard_ID))
summary(body_temp_LMM_bestb)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 |
## date/lizard_ID)
## Data: temp_time_series
##
## REML criterion at convergence: 71392.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8879 -0.3478 -0.0225 0.3300 6.0970
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 1.5384 1.2403
## date (Intercept) 0.3563 0.5969
## Residual 0.9812 0.9906
## Number of obs: 25251, groups: lizard_ID:date, 33; date, 5
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 27.7236 0.4606 11.4180 60.191
## poly(time_min, 2)1 48.0115 1.7051 25212.0897 28.158
## poly(time_min, 2)2 0.6860 1.7158 25212.0585 0.400
## treatmentCooling -2.8157 0.5208 26.3868 -5.406
## treatmentHeating -1.5969 0.5509 27.0887 -2.899
## poly(time_min, 2)1:treatmentCooling -697.0769 2.4336 25213.7698 -286.444
## poly(time_min, 2)2:treatmentCooling 133.9832 2.4104 25212.2645 55.585
## poly(time_min, 2)1:treatmentHeating 484.8482 2.4654 25212.3649 196.660
## poly(time_min, 2)2:treatmentHeating -42.9964 2.4608 25212.2024 -17.473
## Pr(>|t|)
## (Intercept) 1.19e-15 ***
## poly(time_min, 2)1 < 2e-16 ***
## poly(time_min, 2)2 0.68930
## treatmentCooling 1.10e-05 ***
## treatmentHeating 0.00734 **
## poly(time_min, 2)1:treatmentCooling < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001
## ply(tm_,2)2 0.000 -0.077
## tretmntClng -0.584 0.001 0.000
## tretmntHtng -0.560 0.001 0.000 0.480
## ply(_,2)1:C 0.001 -0.701 0.054 0.000 -0.001
## ply(_,2)2:C 0.000 0.055 -0.712 0.000 0.000 -0.005
## ply(_,2)1:H 0.001 -0.692 0.053 -0.001 0.000 0.485 -0.038
## ply(_,2)2:H 0.000 0.054 -0.697 0.000 0.000 -0.038 0.496
## p(_,2)1:H
## ply(tm_,2)1
## ply(tm_,2)2
## tretmntClng
## tretmntHtng
## ply(_,2)1:C
## ply(_,2)2:C
## ply(_,2)1:H
## ply(_,2)2:H -0.021write.csv(broom.mixed::tidy(body_temp_LMM_bestb),
"./Results_Statistics/temp_time_best_model_dorstemp.csv")Rate of Change
We want to know why each lizard had a different rate of change of CEWL during the 15-minute temperature treatment and measurement period.
First, calculate rate of change for each lizard:
change_LM <- lm(data = CEWL_time_series,
CEWL_g_m2h ~
time_min*lizard_ID)
summary(change_LM)##
## Call:
## lm(formula = CEWL_g_m2h ~ time_min * lizard_ID, data = CEWL_time_series)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8164 -0.5799 -0.0359 0.4370 9.2347
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.849726 0.079132 339.303 < 2e-16 ***
## time_min 0.122080 0.008697 14.037 < 2e-16 ***
## lizard_ID402 -0.391770 0.111914 -3.501 0.000465 ***
## lizard_ID404 -4.828469 0.110317 -43.769 < 2e-16 ***
## lizard_ID406 -15.285763 0.118305 -129.207 < 2e-16 ***
## lizard_ID407 -0.098935 0.109533 -0.903 0.366404
## lizard_ID408 -1.653538 0.122070 -13.546 < 2e-16 ***
## lizard_ID409 -6.242781 0.113576 -54.966 < 2e-16 ***
## lizard_ID410 -24.535008 0.111329 -220.384 < 2e-16 ***
## lizard_ID411 -17.712482 0.115498 -153.358 < 2e-16 ***
## lizard_ID412 -16.092335 0.110917 -145.084 < 2e-16 ***
## lizard_ID413 -8.859382 0.113889 -77.789 < 2e-16 ***
## lizard_ID414 -21.472142 0.126188 -170.159 < 2e-16 ***
## lizard_ID415 -17.170442 0.115574 -148.567 < 2e-16 ***
## lizard_ID416 -6.509572 0.111739 -58.257 < 2e-16 ***
## lizard_ID417 -22.062135 0.111838 -197.269 < 2e-16 ***
## lizard_ID418 -18.675743 0.111729 -167.152 < 2e-16 ***
## lizard_ID419 -16.585398 0.112004 -148.079 < 2e-16 ***
## lizard_ID420 -21.060422 0.112011 -188.022 < 2e-16 ***
## lizard_ID421 -20.484787 0.112438 -182.188 < 2e-16 ***
## lizard_ID422 -14.092361 0.114180 -123.422 < 2e-16 ***
## lizard_ID423 -21.460844 0.111367 -192.703 < 2e-16 ***
## lizard_ID424 -19.003214 0.111240 -170.831 < 2e-16 ***
## lizard_ID425 -16.936058 0.111112 -152.424 < 2e-16 ***
## lizard_ID426 -19.507882 0.148579 -131.296 < 2e-16 ***
## lizard_ID427 -14.604353 0.113202 -129.012 < 2e-16 ***
## lizard_ID428 -3.577002 0.120484 -29.689 < 2e-16 ***
## lizard_ID429 -12.050071 0.111385 -108.183 < 2e-16 ***
## lizard_ID430 -13.697643 0.112326 -121.945 < 2e-16 ***
## lizard_ID431 -18.738456 0.111911 -167.441 < 2e-16 ***
## lizard_ID432 -19.967832 0.112583 -177.360 < 2e-16 ***
## lizard_ID433 -9.059766 0.112296 -80.677 < 2e-16 ***
## lizard_ID434 -5.041610 0.112862 -44.671 < 2e-16 ***
## lizard_ID435 -12.863750 0.112050 -114.804 < 2e-16 ***
## lizard_ID436 -5.219401 0.123891 -42.129 < 2e-16 ***
## lizard_ID437 -8.399519 0.112705 -74.527 < 2e-16 ***
## lizard_ID438 -6.781302 0.113911 -59.532 < 2e-16 ***
## lizard_ID439 -11.307001 0.111506 -101.403 < 2e-16 ***
## time_min:lizard_ID402 -0.385888 0.012300 -31.373 < 2e-16 ***
## time_min:lizard_ID404 -0.494876 0.012165 -40.679 < 2e-16 ***
## time_min:lizard_ID406 -0.793039 0.013514 -58.683 < 2e-16 ***
## time_min:lizard_ID407 0.100523 0.012103 8.306 < 2e-16 ***
## time_min:lizard_ID408 -0.561813 0.013056 -43.031 < 2e-16 ***
## time_min:lizard_ID409 -0.622057 0.012433 -50.033 < 2e-16 ***
## time_min:lizard_ID410 -0.020582 0.012247 -1.681 0.092869 .
## time_min:lizard_ID411 -0.463220 0.012420 -37.295 < 2e-16 ***
## time_min:lizard_ID412 0.237357 0.012215 19.432 < 2e-16 ***
## time_min:lizard_ID413 -0.341593 0.012458 -27.419 < 2e-16 ***
## time_min:lizard_ID414 -0.209479 0.017039 -12.294 < 2e-16 ***
## time_min:lizard_ID415 -0.205661 0.012599 -16.324 < 2e-16 ***
## time_min:lizard_ID416 -0.377975 0.012282 -30.775 < 2e-16 ***
## time_min:lizard_ID417 0.291636 0.012292 23.725 < 2e-16 ***
## time_min:lizard_ID418 0.211037 0.012282 17.182 < 2e-16 ***
## time_min:lizard_ID419 0.118183 0.012305 9.604 < 2e-16 ***
## time_min:lizard_ID420 -0.340719 0.012310 -27.679 < 2e-16 ***
## time_min:lizard_ID421 -0.407892 0.012343 -33.047 < 2e-16 ***
## time_min:lizard_ID422 -0.293826 0.012488 -23.529 < 2e-16 ***
## time_min:lizard_ID423 0.567593 0.012255 46.313 < 2e-16 ***
## time_min:lizard_ID424 0.573931 0.012241 46.888 < 2e-16 ***
## time_min:lizard_ID425 -0.183705 0.012230 -15.021 < 2e-16 ***
## time_min:lizard_ID426 -0.820460 0.025602 -32.047 < 2e-16 ***
## time_min:lizard_ID427 -0.013238 0.012403 -1.067 0.285850
## time_min:lizard_ID428 -0.157217 0.012993 -12.100 < 2e-16 ***
## time_min:lizard_ID429 -0.061577 0.012255 -5.025 5.07e-07 ***
## time_min:lizard_ID430 -0.924418 0.012332 -74.963 < 2e-16 ***
## time_min:lizard_ID431 -0.440182 0.012298 -35.793 < 2e-16 ***
## time_min:lizard_ID432 -0.284201 0.012358 -22.998 < 2e-16 ***
## time_min:lizard_ID433 -0.639850 0.012329 -51.897 < 2e-16 ***
## time_min:lizard_ID434 -0.513333 0.012377 -41.473 < 2e-16 ***
## time_min:lizard_ID435 -0.679765 0.012310 -55.219 < 2e-16 ***
## time_min:lizard_ID436 -0.037103 0.015879 -2.337 0.019469 *
## time_min:lizard_ID437 -0.574387 0.012365 -46.452 < 2e-16 ***
## time_min:lizard_ID438 -0.282519 0.012464 -22.667 < 2e-16 ***
## time_min:lizard_ID439 0.538159 0.012266 43.874 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9663 on 28331 degrees of freedom
## Multiple R-squared: 0.9816, Adjusted R-squared: 0.9815
## F-statistic: 2.069e+04 on 73 and 28331 DF, p-value: < 2.2e-16# now get the trend/slope for each individual and make into a variable in df
lizard_coeffs <- data.frame(emtrends(change_LM, "lizard_ID", var = "time_min")) %>%
left_join(read_rds("Data/collection_dat_formatted.RDS"),
by = c("lizard_ID" = "individual_ID")) %>%
dplyr::select(lizard_ID, rate_change = time_min.trend,
SE, df, lower.CL, upper.CL,
date,
treatment,
mass_g, SVL_mm,
chamber_time_elapsed_hr
) %>%
dplyr::mutate(date_factor = as.factor(date))
summary(lizard_coeffs)## lizard_ID rate_change SE df
## 401 : 1 Min. :-0.8023 Min. :0.008416 Min. :28331
## 402 : 1 1st Qu.:-0.3728 1st Qu.:0.008672 1st Qu.:28331
## 404 : 1 Median :-0.1621 Median :0.008739 Median :28331
## 406 : 1 Mean :-0.1074 Mean :0.009524 Mean :28331
## 407 : 1 3rd Qu.: 0.1088 3rd Qu.:0.008920 3rd Qu.:28331
## 408 : 1 Max. : 0.6960 Max. :0.024080 Max. :28331
## (Other):31
## lower.CL upper.CL date treatment
## Min. :-0.81947 Min. :-0.78520 Min. :2021-10-05 Control:12
## 1st Qu.:-0.38947 1st Qu.:-0.35612 1st Qu.:2021-10-11 Cooling:12
## Median :-0.17933 Median :-0.14491 Median :2021-10-12 Heating:13
## Mean :-0.12610 Mean :-0.08877 Mean :2021-10-13
## 3rd Qu.: 0.09151 3rd Qu.: 0.12617 3rd Qu.:2021-10-18
## Max. : 0.67913 Max. : 0.71289 Max. :2021-10-19
##
## mass_g SVL_mm chamber_time_elapsed_hr date_factor
## Min. : 8.70 Min. :60.00 Min. :0.9333 2021-10-05:6
## 1st Qu.:10.00 1st Qu.:64.00 1st Qu.:1.0000 2021-10-11:8
## Median :11.30 Median :66.00 Median :1.0333 2021-10-12:7
## Mean :11.42 Mean :66.24 Mean :1.1005 2021-10-18:8
## 3rd Qu.:12.60 3rd Qu.:69.00 3rd Qu.:1.1500 2021-10-19:8
## Max. :15.30 Max. :75.00 Max. :1.5333
## length(unique(CEWL_time_series$lizard_ID))## [1] 37NOW we can look at how the RATE of CEWL change (CEWL change per MINUTE) is different among lizards.
change_LMM_1 <- lme4::lmer(data = lizard_coeffs,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
mass_g + SVL_mm +
chamber_time_elapsed_hr +
# include date as random effect bc sig diff weather
(1|date_factor))
summary(change_LMM_1)## Linear mixed model fit by REML ['lmerMod']
## Formula: rate_change ~ treatment + mass_g + SVL_mm + chamber_time_elapsed_hr +
## (1 | date_factor)
## Data: lizard_coeffs
##
## REML criterion at convergence: 26.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.97945 -0.60829 0.04974 0.48861 1.63798
##
## Random effects:
## Groups Name Variance Std.Dev.
## date_factor (Intercept) 0.01674 0.1294
## Residual 0.07171 0.2678
## Number of obs: 37, groups: date_factor, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.90245 1.03318 -0.873
## treatmentCooling -0.12485 0.14578 -0.856
## treatmentHeating 0.49210 0.11557 4.258
## mass_g -0.03169 0.03778 -0.839
## SVL_mm 0.01927 0.01781 1.082
## chamber_time_elapsed_hr -0.23232 0.42633 -0.545
##
## Correlation of Fixed Effects:
## (Intr) trtmnC trtmnH mass_g SVL_mm
## tretmntClng 0.194
## tretmntHtng -0.165 0.533
## mass_g 0.396 0.060 -0.084
## SVL_mm -0.880 0.011 0.227 -0.650
## chmbr_tm_l_ -0.346 -0.656 -0.271 -0.157 -0.002car::vif(change_LMM_1)## GVIF Df GVIF^(1/(2*Df))
## treatment 1.934248 2 1.179310
## mass_g 1.825412 1 1.351078
## SVL_mm 1.878111 1 1.370442
## chamber_time_elapsed_hr 1.865277 1 1.365751drop1(change_LMM_1)## Single term deletions
##
## Model:
## rate_change ~ treatment + mass_g + SVL_mm + chamber_time_elapsed_hr +
## (1 | date_factor)
## npar AIC
## <none> 21.929
## treatment 2 42.116
## mass_g 1 20.523
## SVL_mm 1 21.118
## chamber_time_elapsed_hr 1 20.291anova(change_LMM_1)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 2.56693 1.28346 17.8977
## mass_g 1 0.00639 0.00639 0.0891
## SVL_mm 1 0.08371 0.08371 1.1673
## chamber_time_elapsed_hr 1 0.02129 0.02129 0.2969Don’t need to worry about VIF.
drop mass first based on low SS and resulting AIC:
change_LMM_2 <- lme4::lmer(data = lizard_coeffs,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
SVL_mm +
chamber_time_elapsed_hr +
# include date as random effect bc sig diff weather
(1|date_factor))
summary(change_LMM_2)## Linear mixed model fit by REML ['lmerMod']
## Formula: rate_change ~ treatment + SVL_mm + chamber_time_elapsed_hr +
## (1 | date_factor)
## Data: lizard_coeffs
##
## REML criterion at convergence: 22.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8071 -0.6048 0.1496 0.5868 1.7654
##
## Random effects:
## Groups Name Variance Std.Dev.
## date_factor (Intercept) 0.01111 0.1054
## Residual 0.07302 0.2702
## Number of obs: 37, groups: date_factor, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.564781 0.948384 -0.596
## treatmentCooling -0.114887 0.146720 -0.783
## treatmentHeating 0.483926 0.116036 4.170
## SVL_mm 0.009684 0.013537 0.715
## chamber_time_elapsed_hr -0.290802 0.424397 -0.685
##
## Correlation of Fixed Effects:
## (Intr) trtmnC trtmnH SVL_mm
## tretmntClng 0.183
## tretmntHtng -0.139 0.544
## SVL_mm -0.891 0.070 0.224
## chmbr_tm_l_ -0.315 -0.656 -0.291 -0.139drop1(change_LMM_2)## Single term deletions
##
## Model:
## rate_change ~ treatment + SVL_mm + chamber_time_elapsed_hr +
## (1 | date_factor)
## npar AIC
## <none> 20.523
## treatment 2 40.292
## SVL_mm 1 19.118
## chamber_time_elapsed_hr 1 19.059anova(change_LMM_2)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 2.56389 1.28195 17.5558
## SVL_mm 1 0.02860 0.02860 0.3917
## chamber_time_elapsed_hr 1 0.03428 0.03428 0.4695drop SVL:
change_LMM_3 <- lme4::lmer(data = lizard_coeffs,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
chamber_time_elapsed_hr +
# include date as random effect bc sig diff weather
(1|date_factor))
summary(change_LMM_3)## Linear mixed model fit by REML ['lmerMod']
## Formula: rate_change ~ treatment + chamber_time_elapsed_hr + (1 | date_factor)
## Data: lizard_coeffs
##
## REML criterion at convergence: 16.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8677 -0.6281 0.1374 0.5128 1.7759
##
## Random effects:
## Groups Name Variance Std.Dev.
## date_factor (Intercept) 0.01147 0.1071
## Residual 0.07171 0.2678
## Number of obs: 37, groups: date_factor, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.03945 0.42710 0.092
## treatmentCooling -0.12250 0.14505 -0.845
## treatmentHeating 0.46534 0.11207 4.152
## chamber_time_elapsed_hr -0.24832 0.41652 -0.596
##
## Correlation of Fixed Effects:
## (Intr) trtmnC trtmnH
## tretmntClng 0.543
## tretmntHtng 0.136 0.543
## chmbr_tm_l_ -0.977 -0.654 -0.269drop1(change_LMM_3)## Single term deletions
##
## Model:
## rate_change ~ treatment + chamber_time_elapsed_hr + (1 | date_factor)
## npar AIC
## <none> 19.118
## treatment 2 38.317
## chamber_time_elapsed_hr 1 17.510anova(change_LMM_3)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 2.56421 1.28211 17.8795
## chamber_time_elapsed_hr 1 0.02549 0.02549 0.3554drop chamber_time_elapsed_hr:
change_LMM_4 <- lme4::lmer(data = lizard_coeffs,
rate_change ~
treatment +
(1|date_factor))
summary(change_LMM_4)## Linear mixed model fit by REML ['lmerMod']
## Formula: rate_change ~ treatment + (1 | date_factor)
## Data: lizard_coeffs
##
## REML criterion at convergence: 16.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94436 -0.62201 -0.04392 0.39606 1.77788
##
## Random effects:
## Groups Name Variance Std.Dev.
## date_factor (Intercept) 0.01179 0.1086
## Residual 0.07012 0.2648
## Number of obs: 37, groups: date_factor, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.2093 0.0908 -2.305
## treatmentCooling -0.1793 0.1085 -1.653
## treatmentHeating 0.4475 0.1068 4.191
##
## Correlation of Fixed Effects:
## (Intr) trtmnC
## tretmntClng -0.593
## tretmntHtng -0.611 0.504drop1(change_LMM_4)## boundary (singular) fit: see help('isSingular')## Single term deletions
##
## Model:
## rate_change ~ treatment + (1 | date_factor)
## npar AIC
## <none> 17.510
## treatment 2 38.876anova(change_LMM_4)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 2.5645 1.2823 18.286plot(change_LMM_4)
The simplest model is our best/final model, and assumptions look good.
Save with p-values:
change_LMM_best <- lmerTest::lmer(data = lizard_coeffs,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
(1|date_factor))
summary(change_LMM_best)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rate_change ~ treatment + (1 | date_factor)
## Data: lizard_coeffs
##
## REML criterion at convergence: 16.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94436 -0.62201 -0.04392 0.39606 1.77788
##
## Random effects:
## Groups Name Variance Std.Dev.
## date_factor (Intercept) 0.01179 0.1086
## Residual 0.07012 0.2648
## Number of obs: 37, groups: date_factor, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.2093 0.0908 12.1290 -2.305 0.039599 *
## treatmentCooling -0.1793 0.1085 29.9201 -1.653 0.108748
## treatmentHeating 0.4475 0.1068 30.3121 4.191 0.000222 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmnC
## tretmntClng -0.593
## tretmntHtng -0.611 0.504broom.mixed::tidy(change_LMM_best) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/change_CEWL_best_model.csv")
anova(change_LMM_best, type = "3", ddf = "Kenward-Roger")## Type III Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 2.5102 1.2551 2 30.708 17.898 7.034e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Treatment group explained a significant amount of the variation (F(3,30)=18.29, p < 0.0001).
also get pairwise differences for change ~ tmt, which I will use to make the boxplot later:
emmeans <- data.frame(emmeans(change_LMM_best, "treatment"))
emmeans## treatment emmean SE df lower.CL upper.CL
## 1 Control -0.2093260 0.09102635 13.09716 -0.40582833 -0.01282375
## 2 Cooling -0.3886322 0.09189944 12.94683 -0.58725181 -0.19001265
## 3 Heating 0.2381211 0.08845166 11.97139 0.04535039 0.43089180emmean_ps <- data.frame(pairs(emmeans(change_LMM_best, "treatment"))) # p values
emmean_ps## contrast estimate SE df t.ratio p.value
## 1 Control - Cooling 0.1793062 0.1087886 30.33842 1.648208 2.414730e-01
## 2 Control - Heating -0.4474471 0.1074827 30.69471 -4.162968 6.708347e-04
## 3 Cooling - Heating -0.6267533 0.1082740 31.09322 -5.788584 6.511073e-06Figures
COLOR
cool <- c(brewer.pal(11, "RdBu")[c(11)])
heat <- c(brewer.pal(11, "RdBu")[c(2)])
control <- c(brewer.pal(9, "Greys")[c(5)])
cntrls <- c(brewer.pal(9, "Greys")[c(3,6,9,4,7,3,5,8,4,6,9)])cloacal temp ~ time
ggplot(temp_time_series) +
aes(x = time_sec,
y = calibrated_cloacal_temp,
color = treatment) +
geom_point(size = 0.01,
alpha = 0.2) +
geom_smooth(method = "loess",
se = F,
size = 2,
na.rm = TRUE) +
theme_classic() +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time (seconds)" ,
y = "Internal Body Temperature (°C)",
color = "Treatment", se = FALSE) -> ctemp_time## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.ctemp_time## `geom_smooth()` using formula = 'y ~ x'
dorsal temp ~ time
ggplot(temp_time_series) +
aes(x = time_sec,
y = (calibrated_dorsal_temp),
color = treatment) +
geom_point(size = 0.01,
alpha = 0.2) +
geom_smooth(method = "loess",
se = F,
size = 2,
na.rm = TRUE) +
theme_classic() +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x= "Time (seconds)" ,
y= "Surface Body Temperature (°C)",
color ="Treatment") -> dtemp_time
dtemp_time## `geom_smooth()` using formula = 'y ~ x'
CEWL (raw) ~ time
ggplot(CEWL_time_series) +
aes(x= time_sec,
y = CEWL_g_m2h,
color = treatment) +
geom_point(size = 0.01,
alpha = 0.25) +
geom_smooth(method = "loess",
size = 2.5,
se = F) +
theme_classic() +
theme(text = element_text(size = 15)) +
theme(axis.text.x = element_text(size = 10)) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time (seconds)" ,
y = bquote('CEWL (g/'*m^2*'h)'),
color = "Treatment",
se = FALSE) +
scale_x_continuous(limits = c(0, 925),
breaks=seq(0, 950, 100)) -> CEWL_time
CEWL_time## `geom_smooth()` using formula = 'y ~ x'
legend
ggplot(temp_time_series) +
geom_smooth(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
color = treatment),
method = "lm",
formula = y ~ poly(x, 2),
size = 1,
se = F) +
geom_point(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
fill = treatment,
shape = treatment),
size = 4,
alpha = 1) +
scale_x_continuous(limits = c(14, 39),
breaks = c(seq(15, 35, 5)),
labels = c(seq(15, 35, 5))) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = c(0.5,0.5)) +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
scale_fill_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") -> fancy_legend
fancy_legend
LEGEND <- as_ggplot(get_legend(fancy_legend))
# save for using in arranged multi-plots
ggsave(filename = "legend_only.pdf",
plot = LEGEND,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 40, height = 30)CEWL (raw) ~ cloacal temperature
ggplot(temp_time_series) +
geom_point(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
shape = treatment,
color = treatment),
size = 0.01,
alpha = 0.1) +
geom_smooth(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
color = treatment),
method = "lm",
formula = y ~ poly(x, 2),
size = 1.5,
se = F) +
scale_x_continuous(limits = c(12.5, 40),
breaks = c(seq(15, 35, 5)),
labels = c(seq(15, 35, 5))) +
# control group
#annotate(geom = "point", x = 25.4, y = 9.6,
# size = 5, pch = 21, fill = control) +
#annotate(geom = "point", x = 30.14, y = 21.4,
# size = 5, pch = 21, fill = control) +
# annotations for cooling group
annotate(geom = "point", x = 14.3, y = 5,
size = 4, pch = 25, fill = cool) +
annotate("text", x = 12.6, y = 5, label = "T[15]", parse = T,
size = 5, family = "sans", color = cool) +
annotate(geom = "point", x = 35.2, y = 11.8,
size = 4, pch = 25, fill = cool) +
annotate("text", x = 36.7, y = 11.8, label = "T[1]", parse = T,
size = 5, family = "sans", color = cool) +
# annotations for heating group
annotate(geom = "point", x = 15.9, y = 14,
size = 4, pch = 24, fill = heat) +
annotate("text", x = 14.3, y = 14, label = "T[1]", parse = T,
size = 5, family = "sans", color = heat) +
annotate(geom = "point", x = 38.3, y = 25.2,
size = 4, pch = 24, fill = heat) +
annotate("text", x = 40, y = 25.2, label = "T[15]", parse = T,
size = 5, family = "sans", color = heat) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
#legend.text.align = 0,
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
labs(x = "Internal Body Temperature (°C)" ,
y = bquote('CEWL (g '*m^-2*' '*h^-1*')'),
se = FALSE) -> CEWL_ctemp
CEWL_ctemp
# save
ggsave(filename = "CEWL_int_body_temp.pdf",
plot = CEWL_ctemp,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 150, height = 100)CEWL Control
# get control data only
temp_time_series %>%
dplyr::filter(treatment == "Control") %>%
dplyr::filter(complete.cases(calibrated_cloacal_temp, CEWL_g_m2h)) %>%
group_by(lizard_ID) %>%
summarise(max(CEWL_g_m2h)) %>%
nrow()## [1] 11# plot control data only
temp_time_series %>%
dplyr::filter(treatment == "Control") %>%
ggplot() +
aes(x = calibrated_cloacal_temp,
y = CEWL_g_m2h,
color = lizard_ID) +
geom_point(size = 0.1,
alpha = 0.1) +
geom_smooth(method = "lm",
size = 1,
se = F) +
ylim(0,32) +
xlim(25,31) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.text.align = 0,
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
scale_color_manual(values = cntrls) +
labs(x = "Internal Body Temperature (°C)" ,
y = bquote('CEWL (g '*m^-2*' '*h^-1*')'),
color = "Treatment",
se = FALSE) -> CEWL_control
CEWL_control## `geom_smooth()` using formula = 'y ~ x'
# save
ggsave(filename = "CEWL_temp_control_lm.pdf",
plot = CEWL_control,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 80, height = 60)## `geom_smooth()` using formula = 'y ~ x'CEWL (relative) ~ time
# n lizards
CEWL_time_series %>%
dplyr::filter(complete.cases(time_min, relative_CEWL)) %>%
group_by(lizard_ID) %>%
summarise(max(relative_CEWL)) %>%
nrow()## [1] 37# plot
ggplot() +
geom_point(data = CEWL_time_series,
aes(x= time_min,
y = relative_CEWL,
shape = treatment,
color = treatment),
size = 0.01,
alpha = 0.1) +
geom_hline(yintercept = 0, size = 0.5,
linetype = "dashed", color = "black") +
geom_smooth(data = CEWL_time_series,
aes(x = time_min,
y = relative_CEWL,
color = treatment),
method = "lm",
formula = y ~ poly((x), 2),
size = 1.5,
se = F) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none") +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
annotate(geom = "point", x = 15, y = -3.9,
size = 4, pch = 21, fill = control) +
annotate(geom = "point", x = 15, y = -7.5,
size = 4, pch = 25, fill = cool) +
annotate(geom = "point", x = 15, y = 3,
size = 4, pch = 24, fill = heat) +
scale_x_continuous(limits = c(1, 15),
breaks = seq(1, 15, 2)) +
labs(x = "Time (minutes)" ,
y = bquote('CEWL (g '*m^-2*' '*h^-1*')')) -> CEWLr_time
CEWLr_time## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 4 rows containing missing values (`geom_point()`).
# save
ggsave(filename = "CEWL_relative_time.pdf",
plot = CEWLr_time,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 150, height = 100)## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
## Removed 4 rows containing missing values (`geom_point()`).CEWL Rate of Change
nrow(lizard_coeffs)## [1] 37ggplot() +
geom_jitter(data = lizard_coeffs,
aes(x = treatment,
y = rate_change,
shape = treatment,
color = treatment,
fill = treatment),
size = 0.8,
alpha = 0.3,
position = position_jitter(height = 0, width = 0.2)) +
geom_errorbar(data = emmeans,
aes(x = treatment,
y = emmean,
color = treatment,
group = treatment,
ymin = lower.CL,
ymax = upper.CL),
width = .1,
alpha = 0.9) +
geom_point(data = emmeans,
aes(x = treatment,
y = emmean,
shape = treatment,
fill = treatment),
color = "black",
size = 4) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none") +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_fill_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
geom_hline(yintercept = 0, size = 0.5,
linetype = "dashed", color = "black") +
annotate("text", x = 1, y = 0.42, label = "A", size = 3) +
annotate("text", x = 2, y = 0.3, label = "A", size = 3) +
annotate("text", x = 3, y = 0.9, label = "B", size = 3) +
labs(x = "" ,
y = expression(Delta ~ 'CEWL '*min^-1*''),
color = "Treatment") -> rate_change_CEWL_fig
rate_change_CEWL_fig
# save
ggsave(filename = "CEWL_change_min.pdf",
plot = rate_change_CEWL_fig,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 80, height = 60)Experimental Design Skematic
exp_design <- read.csv("./Data/exp design.csv",
header = TRUE,
fileEncoding="UTF-8-BOM")ggplot(exp_design) +
aes(x = time_min,
y = temp_C,
color = treatment) +
geom_line(size = 1) +
geom_vline(xintercept = 60,
alpha = 0.5,
size = 0.5,
linetype = "dashed", color = "black") +
geom_text(aes(30, 37, family = "sans"), label = "a", color = "black", size = 3) +
geom_vline(xintercept = 65,
alpha = 0.5,
size = 0.5,
linetype = "dashed", color = "black") +
geom_text(aes(62.5, 37, family = "sans"), label = "b", color = "black", size = 3) +
geom_vline(xintercept = 80,
alpha = 0.5,
size = 0.5,
linetype = "dashed", color = "black") +
geom_text(aes(72.5, 37, family = "sans"), label = "c", color = "black", size = 3) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
axis.text.x = element_blank(),
#axis.line.x = element_blank(),
axis.ticks.x = element_blank()) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time" ,
y = "Body Temperature (°C)",
color = "Treatment") +
annotate(geom = "point", x = 80, y = 25,
size = 3, pch = 21, fill = control) +
annotate(geom = "point", x = 80, y = 20,
size = 3, pch = 25, fill = cool) +
annotate(geom = "point", x = 80, y = 30,
size = 3, pch = 24, fill = heat) +
scale_y_continuous(limits = c(14, 37),
breaks=seq(15, 35, 5)) -> methods_schmatic_fig
methods_schmatic_fig
# save
ggsave(filename = "methods_schematic.pdf",
plot = methods_schmatic_fig,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 80, height = 60)Save Grouped Figures
Relative CEWL & Rate of Change
multi_fig_relative <- ggarrange(CEWLr_time,
ggarrange(rate_change_CEWL_fig, LEGEND,
nrow = 1, ncol = 2,
widths = c(2,1),
align = "hv"
),
nrow = 2, ncol = 1,
heights = c(2,1),
labels = c("A", "B"))## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 4 rows containing missing values (`geom_point()`).multi_fig_relative
ggsave(filename = "CEWL_time_multi_fig.pdf",
plot = multi_fig_relative,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 120, height = 160) # full width = 180Raw CEWL ~ Cloacal temp & control zoom-in
multi_fig_CEWL_by_temp <- ggarrange(CEWL_ctemp,
ggarrange(CEWL_control, LEGEND,
nrow = 1, ncol = 2,
widths = c(2, 1),
align = "hv"
),
nrow = 2, ncol = 1,
heights = c(2,1),
labels = c("A", "B"))## `geom_smooth()` using formula = 'y ~ x'multi_fig_CEWL_by_temp
ggsave(filename = "CEWL_temp_multi_fig.pdf",
plot = multi_fig_CEWL_by_temp,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 120, height = 160)